# SI Figure 1
rm(list=ls())
gc()

quartz(width = 10, height = 5, type = "pdf", "SI Figure damages", file = "SI_Figure_damages.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,1,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(3,3), heights=c(3,3))

# Panel 1: Reflecting temperature distribution through damage function
# Damage distributions
t <- seq(0,10,0.01)
DN <- (1-(1/(1+0.0023*t^2)))
DW <- (1-(1/(1+(t/20.46)^2 + (t/6.081)^6)))

plot(t,DN,type="n",xlim=c(0,10),ylim=c(0,1),xaxs="i",yaxs="i",bty="L",yaxt="n",xaxt="n",xlab="",ylab="")
lines(t,DN)
lines(t,DW)
axis(side=1,at=seq(0,10,2))
axis(side=2,at=seq(0,1,.2),las=2)

lines(x=c(8,8), y=c(0,DW[t==8]), lty=2)
lines(x=c(0,8), y=c(DN[t==8],DN[t==8]), lty=2)
lines(x=c(0,8), y=c(DW[t==8],DW[t==8]), lty=2)

round(DW[t==8]*100,0)
round(DW[t==4]*100,0)
text(x=0, y = DW[t==8]+.03, labels=expression("Weitzman's function:"),pos=4,cex=.8)
text(x=0, y = DW[t==8]-.05, labels=expression("84% of output lost when " * Delta * "T=8,"),pos=4,cex=.8)
text(x=0, y = DW[t==8]-.10, labels=expression("8x compared with " * Delta * "T=4."),pos=4,cex=.8)

round(DN[t==8]*100,0)
text(x=0, y = DN[t==8]+.03, labels=expression("Nordhaus' function:"),pos=4,cex=.8)
text(x=0, y = DN[t==8]-.05, labels=expression("13% of output lost when " * Delta * "T=8."),pos=4,cex=.8)

mtext(side=1, expression("Global mean temperature anomaly, " * Delta * "T (˚C)"), outer=T, line=2, cex=.9, at=0.5)
mtext(side=3, expression("Share of output"), outer=T, line=-3, cex=.9, at=0.1)
mtext(side=3, expression("lost, 1-D(" * Delta * "T)"), outer=T, line=-4, cex=.9, at=0.1)


# Panel 2: Reflecting temperature distribution through expected utility function
# Economic assumptions
eta = 1.45 # Marginal utility of consumption
gdp = 80000 # Aggregate consumption (billion USD)
pop = 8 # Population (billion people)
gdp.pc <- (gdp/pop)/1000 # Undistured per capita consumption (thousand USD).

# Utility scalars
UofC <- ((gdp.pc)^(1-eta))/(1-eta) # Utility of undisturbed consumption
scalar <- UofC-(((gdp.pc-0.001)^(1-eta))/(1-eta)) # Utility value of the marginal dollar

## Utility loss function
EUN <- (UofC - ((gdp.pc*(1/(1+0.0023*t^2)))^(1-eta))/(1-eta))/scalar
EUW <- (UofC - ((gdp.pc*(1/(1+(t/20.46)^2 + (t/6.081)^6)))^(1-eta))/(1-eta))/scalar

plot(t,t,type="n",xlim=c(0,10),ylim=c(0,max(EUW)),xaxs="i",yaxs="i",bty="L",yaxt="n",xlab="",ylab="")
lines(t,EUN)
lines(t,EUW)

lines(x=c(8,8), y=c(0,EUW[t==8]), lty=2)
lines(x=c(0,8), y=c(EUN[t==8],EUN[t==8]), lty=2)
lines(x=c(0,8), y=c(EUW[t==8],EUW[t==8]), lty=2)

EUW[t==8]/EUW[t==4]
text(x=0, y = EUW[t==8]+9000, labels=expression("Weitzman's function:"),pos=4,cex=.8)
text(x=0, y = EUW[t==8]+5000, labels=expression("Utility losses when " * Delta * "T=8 are"),pos=4,cex=.8)
text(x=0, y = EUW[t==8]+1500, labels=expression("25x compared with " * Delta * "T=4."),pos=4,cex=.8)

EUN[t==8]/EUN[t==4]
text(x=0, y = EUN[t==8]+9000, labels=expression("Nordhaus' function:"),pos=4,cex=.8)
text(x=0, y = EUN[t==8]+5000, labels=expression("Utility losses when " * Delta * "T=8 are"),pos=4,cex=.8)
text(x=0, y = EUN[t==8]+1500, labels=expression("4x compared with " * Delta * "T=4."),pos=4,cex=.8)

mtext(side=3, "Utility losses", outer=T, line=-3, cex=.9, at=0.6)
mtext(side=3, expression("u(D(0)) - u(D(" * Delta * "T))"), outer=T, line=-4, cex=.9, at=0.6)

axis(side=2,at=c(0,max(EUW)),label=c("",""),las=2)


dev.off()